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I . REPORT  SUMMARY 

This  is  the  final  quarterly  technical  report  on  a one-year  contract 
to  extend  and  improve  an  existing  empirical  model  for  worldwide  behavior 
of  ionosphei  i.cally  imposed  radio-wave  scintillation.  The  objectives  of 
the  project  were  (1)  to  improve  the  accuracy  of  model-based  calculations 
of  the  intersity-scintillation  index  and  (2)  to  develop  a capability  for 
full  description  (from  the  point  of  view  of  engineering  a /plications)  of 
the  first-order,  complex-signal  statistics  that  characterize  the  trans- 
ionospheric  radio  communication  channel.  A follow-on  project  has  been 
initiated  to  extend  the  channel  model  to  include  second-order  signal 
statistics  in  the  temporal,  spatial,  and  spectral  domains.  The  first 
priority  in  the  follow-on  work  'ill  be  to  calculate  the  fluctuation 
spectra  of  relevant  signal  parameters.  In  addition,  it  is  intended  to 
extend  validity  of  the  model  into  the  multiple-scatter  regime# 


The  stated  objectives  of  the  just-completed  project  were  accomplished 
by  starting  with  the  scintillation  mod  1 developed  by  Fremouw  and  Rino 
(1973),  adding  a geomagnetic-activity  dependence  for  the  behavior  of 
scintillation-producing  ionospheric  irregularities,  and  substituting  a 
more  accurate  and  general  scattering  theory  for  that  of  Briggs  and  Parkin 
(1963)  used  in  the  earlier  model.  The  results  of  these  efforts  have 
been  implemented  in  two  computer  programs.  The  first,  entitled  RFMOD, 
contains  the  main  elements  of  the  scattering  theory,  the  morphological 
model  for  ionospheric  irregularity  strength  and  other  parameters,  and 
subroutines  for  calculating  the  geometry  and  other  relevant  quantities 
from  user-specified  inputs.  The  main  outputs  of  RFMOD  are  the  intensity- 
scintillation  index  S^,  which  is  the  fractional  rms  fluctuation  in  re- 
ceived signal  intensity  (Briggs  and  Parkin,  1963),  and  an  analogous  phase 
scintillation  index  (labled  PRMS) , which  is  the  rms  fluctuation  of 

phase . 
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Three  related  first-order,  signal-statistical  parameters  also  are 
output  from  RFHOD , to  be  used  as  inputs  to  the  second  program,  called 
DIST.  The  latter  program  permits  calculation  of  the  first-order  dis- 
tribution functions  of  amplitude  and  phase  for  selected  values  of 
and/or  $ rmS . There  are  three  types  of  output  from  DIST:  values  of  the 
probability  density  functions  (pdf)  for  amplitude  and  for  phase,  and 
values  of  the  cumulative  probability  function  (cpf)  for  intensity  (or, 
equivalently,  for  amplitude  since  the  argument  is  output  in  dB  relative 
to  the  undisturbed  level).  The  cpf  output  constitutes  estimation  of 
fade  margin  required  to  combat  amplitude  scintillation. 

RFMOD  and  DIST  represent  interim  versions  of  evolving  transiono- 
spheric  communication-channel  codes.  As  stated  above,  capabilities 
will  be  added  for  calculating  second-order  statistics.  More  relevant 
in  the  near  term  is  that  there  has  been  no  direct  testing  of  calculations 
relating  to  phase  scintillation.  As  was  recognized  at  the  outset  of  the 
project,  data  are  not  available  at  present  that  are  directly  relevant 
to  testing  phase  statistics  unambiguously.  Useful  data  are  expected 
from  forthcoming  observations  of  mul ti frequency  coherent  satellite 
beacons,  and  such  tests  will  be  made.  The  forthcoming  data  may  also 
result  in  model  changes  that  will  further  alter  the  calculated  results 
for  amplitude  scintillation,  because  testing  to  date  has  still  necessarily 
relied  on  some  observations  reported  in  terms  of  a partly  subjective 
scintillation  index.  In  addition,  while  the  assumption  of  weak,  single 
scatter  contained  in  the  Briggs-and-Parkin  theory  has  been  partially 
removed,  a more  complete  correction  for  multiple  scatter  effects  has  yet 
to  be  incorporated. 

This  report  represents  a summary  of  the  work  performed  during  the 
past  year  and  is  inteded  as  a status  report  on  the  evolving  SRI/ARPA 
channel  model.  Section  II  describes  the  scattering  theory  employed  and 

I > 
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compares  some  of  its  salient  points  with  those  of  other  theories.  Ihat 
section  includes  some  material  reported  in  earlier  quarterly  reports, 
for  the  sake  of  presenting  a coherent  discussion  under  one  cover,  but 
there  is  no  attempt  to  describe  the  theoretical  development  in  chrono- 
logical detail.  Section  III  contains  a brief  review  of  the  known  mor- 
phology of  scintillation-producing  ionospheric  irregularities  as  modeled 
in  RFMOD . A functional  description  of  the  programs  RFMOD  and  DIST  is 
presented  in  Section  IV,  along  with  instructions  for  their  use.  Finally, 
in  Section  V,  we  present  our  conclusions  on  the  utility  and  limitations 
of  the  channel  model  in  its  present  form  and  some  prospects  for  further 
improvement . 


II  SCATTERING  THEORY  EMPLOYED 

The  basic  approach  to  modeling  the  transionospheric  communication 
channel  is  to  treat  the  channel  as  a time -and -space -varying  linear  filter 
linking  the  transmitting  and  receiving  antennas  of  a communication  sys- 
tem. Slowly  varying  effects  (on  a time  scale  of  many  minutes  or  hours), 
such  as  group  delay,  polarization  change,  and  dispersion  of  signals 
propagating  through  the  smooth  ionosphere,  can  be  described  by  well- 
known  deterministic  methods  and  are  not  included  in  this  work.  What  is 
treated  are  the  relatively  rapid  and  random  variations  (on  a time  scale 
of  fractions  of  a second  to  a minute  or  so)  in  signal  parameters  that 

arise  due  to  scattering  by  irregularities  in  ionospheric  electron  density  I 

and  that  are  referred  to  as  scintillations.  Because  the  channel  can  be 
treated  as  a linear  filter,  standard  Fourier  techniques  can  be  used  to 
apply  the  results  to  modulated  signals  if  a sufficiently  general  descrip- 
tion is  provided  of  the  signal  statistics  resulting  when  a continuous 
(CW)  is  passed  through  the  channel  (Fremouw,  1969).  In  this  section, 
we  discuss  the  first-order  statistics  of  this  basic  signal. 
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Any  cW  signal  can  be  represented  by  a vector  on  the  complex  plane, 
such  as  the  phasor  E in  Figure  1.  In  this  figure,  E = ( E^_  + iE^  J is  meant 
to  represent  the  randomly  time-varying  complex  amplitude  of  the  signal 
output  from  an  antenna  terminal.  In  describing  the  first-order,  complex- 
signal  statistics  of  E,  we  are  assessing  the  probability  that  the  tip  of 
the  phasor  lies  within  an  elemental  area  on  the  complex  plain.  To  begin, 
we  write  E as  the  resultant  of  its  long-term  mean  E and  a zero-mean, 
randomly  varying  component  E . We  choose  to  reference  all  phases  to 
that  of  (e),  which  is  the  same  as  the  phase  that  E would  have  in  the 


LA  341 6-12 

FIGURE  1 PHASOR  DIAGRAM  FOR  CONSIDERATION  OF  FIRST-ORDER 
SIGNAL  STATISTICS,  SHOWING  ELLIPTICAL  CONTOUR  OF 
EQUAL  PROBABILITY  PREDICTED  BY  THE  GAUSSIAN- 
STATISTICS  HYPOTHESIS 


absence  of  channel -imposed  perturbations  because  we  have  excluded  slow 
channel  variations  from  the  model.  (Thus  (E  =0,  and  E = i = 

y x 

<E>.)  The  magnitude  of  (E>  (=  T])  is  less  than  the  amplitude  of  the  un- 
disturbed signal,  however,  because  some  of  the  transmitted  energy  has  been 
scattered  into  the  randomly  varying  component  E . 
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If  it  should  be  that  the  quadrature  components,  E and  E , are 

x y 

jointly  Gaussian  random  variables,  the  utility  of  the  three  variances 

defined  above--or  alternatively  of  and  B^'-would  be  very  great.  Indeed, 

in  this  instance  the  first-order  statistics  of  the  received  signal  E are 

totally  defined  by  either  of  these  sets  of  three  parameters  (i.e.,  the 

three  real  variances  or  the  real  R^  and  the  real  and  imaginary  par^s  of 

B ) . Most  scintillation  observations  have  been  only  in  terms  of  the 
o 

magnitude  of  E (i.e.,  amplitude  or  intensity  scintillation)  because  a 
sufficiently  stable  phase  reference  for  measuring  its  complex  value  has 


seldom  been  available  under  scintillation  conditions.  Thus,  there  has 
been  no  direct  determination  of  the  complex-signal  statistics.  (This 


data  deficiency  should  be  remedied  within  the  next  year  by  means  ot 
careful  measurements  of  signals  from  the  VHF-UHF  coherent  beacons  on 
ATS-6  and  the  Transit  satellites  and,  under  a wider  range  of  ionospheric 
conditions,  by  means  of  observations  of  the  forthcoming  DNA -002  multi- 
frequency beacon.) 

The  Gaussian-statistics  hypothesis  (i.e.,  that  the  quadrature  com 
ponents  are  jointly  Gaussian  variates)  has  been  tested  indirectly,  however, 
by  means  of  amplitude-scintillation  data.  Rino , Livingston,  and  Whitney 
(1975)  assessed  this  hypothesis  and  the  competing  one  of  log-normal 
statistics  (according  to  which  the  phase  and  the  logarithm  of  amplitude 
are  jointly  Gaussian  variates)  by  performing  chi-square,  goodness-of-f it 
tests  of  amplitude  histograms  to  amplitude  pdf's  calculated  from  the  two 
hypotheses.  Good  fits  were  found  in  general  for  both  hypotheses,  which 
tend  toward  identical  pdf's  as  the  scintillation  index  deceases,  but 
the  Gaussian  hypothesis  yielded  better  fits  to  11  out  of  the  12  data 
sets  tested.  This  result  is  fortunate  because  it  permits  presentation 
of  scintillation  calculations  in  a form  that  will  be  rather  familiar  to 
communications  engineers.  In  the  remainder  of  this  report,  we  proceed 
on  the  basis  of  the  Gaussian  signal-statistics  hypothesis. 

For  Gaussian  signal  statistics,  a contour  of  equal  probability  for 
the  tip  of  the  E-phascr  on  the  complex  plane  is  an  ellipse  (Beckman  and 
Spizzichino,  1963),  such  as  the  one  in  Figure  1.  The  analytical  task  of 
first-order,  signal-statistical  channel  modeling,  then,  is  to  relate  the 
parameters  of  such  an  ellipse  to  parameters  of  the  channel  (scattering 
irregularities  in  the  ionosphere  in  the  present  case).  The  probability 
ellipse  is  characterized  by  its  size,  eccentricity,  and  orientation  on  the 
complex  plane.  Clearly,  Rq  (being  the  fraction  of  flux  scattered)  dic- 
tates the  size  of  the  ellipse.  The  eccentricity  is  dictated  jointly  by 
R and  the  magnitude  of  Bq,  and  the  orientation  is  controlled  by  the  phase 
angle  of  B . Spec: fic illy , the  axial  ratio  of  the  probability  ellipse  is 

given  by 


1/2 


r 


(3) 


and  its  orientation  angle,  6 (shown  on  Figure  1),  is  given  by 


6 


(4) 


It  will  be  recognized  from  Eqs.  (1)  through  (4)  and  the  accompanying 

discussion  that  there  are  several  sets  of  three  parameters  that  can  be 

used  to  describe  the  first-order,  complex-signal  statistics.  For  reasons 

of  convenience,  which  will  become  clear  in  the  ensuing  discussion,  we 

choose  to  work  with  the  parameters  R , |b  |/R  , and  6.  As  an  example  of 

o o o 

the  relationship  of  these  three  parameters  to  a quantity  of  direct 

applications  interest,  we  have  plotted  the  intensity-scintillation  index 

S in  Figure  2 from  the  following  formula  derived  by  Nakagami  (1960) : 

4 


o 


We  note  that  in  Briggs'  and  Parkin's  theory,  is  calculated 
approximately  as 


S 


2 

4 


o 


(6) 
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FIGURE  2 SCINTILL  ATION  INDEX  S4  AS  A FULCTION  OF  ELONGATION  AND 
ORiENTATION  OF  THE  EQUAL-PROBABILITY  ELLIPSE,  FOR  FOUR 
VALUES  OF  IONOSPHERIC  SCATTERING  COEFFICIENT,  Rq 


The  following  three  relationships  arise  from  Eqs . (1)  and  (2): 


„ R |B  | 

2 o . I o I ,„cs 

o = — 1 - cos  (2o) 

x 2 R 

o 


r r |B  i 1 

0^  = — I 1 + ° cos  (26) 


1 


2 

a 

xy 


o 


(26) 


(9) 


Thus,  Eq.  (6)  contains  o .ly  the  contribution  of  the  in-phase  scattered 
component  S to  S . It  arises  from  equating  fluctuations  in  E^  to 
fluctuations  in  the  real  amplitude  A and  is  accurate  orJ  y for  very  weak 
scatter  [i.e.,  when  the  higher-order  terms  in  Eq.  (5)  are  negligible]. 

A significant  improvement  in  the  present  scintillation  model  is  retention 
of  all  terms  in  Eq.  (5)  , which  accounts  for  the  contribution  of  fluctua- 
tions in  the  phase-quadrature  component  E^(=S^)  to  amplitude  scintilla- 
tion (as  well  as  for  che  influence  of  the  correlation  between  and  S ) . 

Other  differences  between  the  old  scattering  theory  and  the  new 

arise  in  the  manner  of  calculating  R and  B . The  basic  approach  used 

o o 

in  the  present  work  is  that  described  by  Rino  and  Fremouw  (1973b),  in 
which  the  scattering  geometry  is  treated  quite  rigorously.  Briggs  and 
Parkin  simplified  the  geometry  by  assuming  that  the  ionospheric  scattering 
layer  effectively  can  be  tilted  from  the  horizontal  so  that  its  boundaries 
are  always  normal  to  the  incident  radio-wave  propagation  vector.  This 
approximation  breaks  down  for  large  zenith  angles  but  it  avoids  a great 
deal  of  complication  in  derivation  of  the  scattering  equations.  Some  of 
this  complication  was  avoided  by  Rino  and  Fremouw  (1973b)  by  assuming 
zero  correlation  in  electron-density  fluctuation  along  the  vertical,  which 
has  negligible  effect  on  the  scattering  equations  for  nearly  vertical 
incidence.  Again,  however,  the  effect  of  the  assumption  becomes  un- 
acceptable as  the  incidence  angle  increases. 

During  the  last  quarter  of  this  contract,  the  theory  was  recast  to 
include  the  scattering  effect  of  structure  in  all  three  dimensions  of 
the  ionosphere.  An  outline  of  the  derivation  is  given  here,  together 
with  the  resulting  equations  for  describing  first-order  signal  statistics, 
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which  are  valid  over  a greater  range  of  incidence  angles  than  their 
counterparts  in  Quarterly  Technical  Reports  1 and  2. 

The  scattering  coefficient  represents  a special  case  of  the 
spatial  autocorrelation  function  of  the  scattered  signal,  which  is 
defined  as 

R(P)  = <E  (r)  e"  (r  + P) ) ; (10) 

s s 

namely,  the  case  in  which  the  spatial  lag  parameter  p is  zero.  A similar 
definition,  calculated  without  taking  the  complex  conjugate,  can  be  made 
of  a quantity  B(p) , which  reduces  in  the  same  special  case  (i.e., 
zero  lag) . The  approach  to  evaluating  R(p)  and  B(p)  is  to  express  Eg 
as  the  scattered  field  calculated  by  means  of  the  first  Born  approxi- 
mation, with  the  scattering  medium  represented  by  a spatial  spectrum  of 
ionospheric  structure.  The  field  is  calculated  by  integrating  the 
contribution  from  different  heights  within  the  scattering  layer,  con- 
sidering the  diffraction  effects  arising  in  propagation  to  the  observing 
plane  but  ignoring  effects  of  multiple  scatter.  When  the  complex  field 
is  multiplied  by  its  complex  conjugate,  the  effects  cancel. 

The  result  is  that  R(o)  can  be  expressed  as  the  two-dimensional 
Fourier  transform  of  a height-smeared  version  of  the  ionosphere's  spatial 
spectrum  and  B(p)  can  be  expressed  as  the  transform  of  that  spectrum 
times  a propagation  factor.  Setting  the  lag  parameter  to  zero  is  tanta- 
mount to  integrating  over  the  spectra  in  the  observing  plane  (as  opposed 
to  calculating  the  Fourier  transforms),  with  being  an  integral  over 
the  scattered  field's  angular  power  spectrum  and  B^  being  an  integral 
over  a complex  spectrum  that  incorporates  diffraction  effects. 

krobably  more  important  for  first-order  statistics  than  the  approacli 
to  deriving  R and  B is  the  form  of  the  spatial  spectrum  used  in 


10 


their  calculation. 


In  evaluating  R and  the  second  factor  in  Eq.  (6), 

o 

Briggs  and  Parkin  assumed  a Gaussian-shaped  spatial  spectrum  for  the 
ionosphere.  While  convenient  because  of  its  Fourier-transformation 
characteristics,  a Gaussian  spectrum  does  not  appear  to  represent 
ionospheric  structure  very  realistically  (Rufenach,  1^72;  Dyson  et  al . , 
1974)  . 

A better  description  is  afforded  by  a power  law  with  a low-frequency 

cutoff.  When  normalized  so  that  its  three-dimensional  integral  is  unity 

2, 

times  the  spatial  variance  <(AN)  ) of  electron  density,  such  a spectrum 

3 ^ 

[units  of  (electrons/m  ) per  unit  three-dimensional  spatial  frequency, 

—4 

*3  may  be  written  as 


3/2  r(v  + 1/2) 
T(v  - 1) 


c^a  ((AN)  ) 

- 2 v + 1/2 
1 + [aS(*)j 


where  v is  the  spectral  index  defining  the  power  law,  a is  scale-size 
parameter  controlling  the  low-frequency  cutoff,  a is  the  ratio  of 
irregularity  scale  size  along  the  geomagnetic  field  to  that  across  the 
field,  and  the  three  dimensional  shape  and  orientation  >f  the  irregulari- 
ties (Budden,  1965),  are  described  by 


S2(K)  = A(f)  K2  + K2  + BU)*2  - 2C(^)K  K 
x y z x z 


C(t)  = (1  - a ) cos  4 sin  i|i 


(15) 


and  4 = the  magnetic  dip. 

The  coordinate  system  used  is  such  that  K is  directed  to  the  geomagnetic 
north  K to  the  geomagnetic  east,  and  it  downward. 

y z 

Using  the  above  spectrum,  one  obtains  the  following  expression  for 


R : 
o 


r-  2 . 2 T(V  - 1/2)  a a sec __  6 , //Ax„2N 

R = 2yn  r X “fT 77^  ~ L((AN) 

o e I (v  - 1)  A 'c 7 


(16) 


where 


and 


r = classical  electron  radius, 
e 

X = observing  wavelength, 

L = thickness  of  scattering  layer, 

A'  = 1/2[d'  + D], 

C'  = 1/2[d'  - D], 


2 

D ' = 1 + B ( ^ + 90°)  + 2C(>ii)  tan  0 cos  0 + B(i{i)  tan  0, 


D = ./[O'  - 2(1  + sin^  0)  p-  + (B  )^  , 

2 

b'  = B(4')  tan  9 sin  20  + 2C(4)  tan  0 sin  0, 

» 

0 = incidence  angle  on  scattering  layer, 

0 = propagation  azimuth  relative  to  geomagnetic  north  at 

height  of  scattering  layer. 


Eq.  (16)  was  obtained  by  integrating  analytically  over  the  angular 

power  spectrum  of  the  scattered  wave  arriving  at  the  observing  plane, 

as  described  above.  The  integral  over  the  complex  spectrum  associated 

with  B cannot  be  performed  analytically.  However,  the  real  and  imaginary 
o 

integrals  involved  have  been  solved  numerically  for  a large  range  of 
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geometrical  situations  described  by  the  parameters  a,  and  Study 

of  the  numerical  results  has  led  to  approximations  suitable  for  implemen- 
tation in  an  operating  computer  code  such  as  RFMOD . 

Simplifications  have  been  found  by  converting  the  real  and  imaginary 
integrals  to  the  magnitude  and  phase  angle  of  Bq  and  then  to  the  quantities 
|B  l/R  and  6 [the  latter  via  Eq.  (4)J.  The  simplification  arises  from 
an  empirically  noted  approximate  relationship  between  the  latter  two 
quantities,  as  described  in  detail  in  Quarterly  Technical  Reports  1 and  2. 
It  permits  approximate  description  of  the  complete  first-order,  complex- 
signal  statistics  from  solution  of  the  real  integral  alone.  What  is 
actually  calculated  is  C*/Rq,  which  is  related  to  ReilM  by  Eq.  (7). 

Two  approximations  are  involved  in  formulating  this  part  of  the 
problem.  The  first-neglect  of  a diffraction  effect  related  to  layer 
thickness --appears  numerically  to  be  valid  essentially  anywhere  outside 
the  scattering  layer.  The  second  holds  so  long  as  the  spectrum  of 
irregularities  contributing  to  E^  fluctuations  is  cut  off  at  the  low- 
frequency  end  by  diffraction  effects  rather  than  by  an  inherent  cutoff 
of  the  ionospheric  spatial  spectrum.  Both  in-situ  (Dyson  et  al.,  1974) 
and  scintillation  observations  (Rufenach,  1972;  Rino,  Livingston,  and 
Whitney,  1975)  support  this  view  for  all  frequencies  of  interest.  We 
have,  therefore,  proceeded  on  the  basis  of  this  "near-zone"  approximation, 
but  we  have  incorporated  a test  in  RIT10D  to  safeguard  against  violation 
of  the  assumption  underlying  it.  On  the  basis  of  these  two  approximations 
2 

one  may  write  O /R  as  follows: 

7 x o 

2 , — — / \ v - i / 2 

\ = 2(v  - 1/2)  VA  AZ  sec  6\ 

R0  (c')V\4tt  a2  / 
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(17) 


°°  2tt 

ri 


sin2  [q2  7(9 ,$') 


r i 2 2 . 2 nv  + 1/2 

L(P  ) cos  9 + sm  9-1 


d9  dq 


where  \A.Z  = Fresnel-zone  radius  at  scattering  layer 

2 2 

7(9,0')  = [l  + tan  6 cos  (0'  - 9)3 

P'  = v/a'/c' 

< * 1 -1 
and  0 = 0 - ~ tan 


D'  - 2(1  + sin  0) 


(18) 


(19) 


The  main  diffraction  effect  dictating  the  first-order  statistics  of 
a scattered  signal  is  controlled  by  the  Fresnel-distance  parameter 
appearing  in  parentheses  in  front  of  the  integral  in  Eq.  (17).  The 
integral  itself  contains  only  geometric  factors.  Although  the  integral 
cannot  be  solved  analytically,  it  has  been  possible  to  find  a useful 
algebraic  approximation  to  it.  The  approach  was  to  numerically  evaluate 
the  integral,  which  was  named  F(P',9,  0,),  for  selected  points  in  a 
three-space  of  P',  9,  and  0',  and  then  to  fit  judiciously  chosen 
algebraic  functions  to  the  points.  The  resulting  functions  then  were 
incorporated  into  RFMOD.  The  procedure  was  as  follows. 

First,  the  integral  was  solved  numerically  for  p = 1,  2,  5,  10,  and 
20,  with  9 ranging  from  0 through  80°  (which  exceeds  the  maximum  in- 
cidence angle  of  a trans ionospheric  communication  link  on  the  F layer, 
by  virture  of  the  ionosphere's  curvature),  and  with  0 stepped  between 
0°  and  90°  (since  the  0*  dependence  is  symmetrical  about  90  ).  After 
the  results  were  plotted,  inspection  of  the  curves  suggested  a behavior 
of  the  following  form: 

F(P\  9,  0')  « ~ a ( P r , 9)  + b(P',  9)  + (b  - a)  cos  20',  (20) 
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where  a = the  value  of  the  integral  when  1 - TT/2,  a \d  b t ie  value 

whe  n 1 ' = 0 . 

The  values  of  a and  b were  then  plotted  as  functions  of  0 for  given 
values  of  3'.  Again  the  plots  suggested  analytic  fits  of  the  form 


a(p',e)  = c(g  ')  sec  ft 


m(B  ') 


b(P  ' ,6)  = c(P  ')  sec  0 


n(P') 


(21) 

(22) 


where  c(P ')  = the  value  of  the  function  when  0 = 0 =0  and  where  the 

exponents  m(P ')  and  n(P  ')  were  chosen  by  forcing  Eqs . (21)  and  (22)  to 
match  the  numerical  value  of  the  integral  when  0 = 70  , which  is  about 
its  largest  value  for  penetration  of  the  F layer.  The  values  of  c,  m, 
and  n were  then  employed  in  a one-dimensional,  least-squares  fitting 
routine  to  obtain  best-fit  polynomial  exaansio.."  of  c(P  ) , m(P  ) , and 
n(g  ')  . 

The  results  of  the  above  procedure  were  combined  into  an  algebraic 
approximation  to  the  integral  in  Eq.  (17),  and  the  approximation  was 
then  used  to  calculate  curves  of  F(P',  6,  1 ')  over  the  ranges  of  0 and 
l'  and  for  the  values  of  originally  employed  in  the  numerical  in- 
tegration. Finally,  similar  curves  were  calculated  for  other  (both 
intermediate  and  more  extreme)  values  of  P , and  numerical  integrations 
were  performed  for  these  new  cases.  The  entire  collections  of  numerical- 
integration  and  algebraic-approximation  results  were  then  compared  to 
determine  the  efficacy  of  the  approximate  expression.  An  example  of  the 
comparison  is  presented  in  Figure  3 for  a value  of  P not  involved  in 
development  of  the  approximation;  the  smooth  curves  represent  the 
algebraic  fit  and  the  individual  points  are  results  of  the  numerical 
integration.  The  fit  is  considered  quite  adequate  for  application  to 
scintillation  modeling.  Thus,  RFMOD  contains  the  following  expression: 


I 


solution  to  the  imaginary  part  of  the  Bo  integral  has  been  found,  and 

numerical  integration  each  time  RFKOD  is  run  seems  unacceptably  inefficient, 

Fortunately,  as  mentioned  above,  synthesis  of  the  numerical  integration 

results  that  have  been  obtained  has  pointed  to  another  simple  procedure. 

Specifically,  when  |b  |/R  and  6 were  formed  from  the  real  and  imaginary 

integrals  for  a number  of  geometrical  situations,  it  was  found  that  tie 

relationship  between  |b  | /R  and  6 was  not  very  dependent  on  geometry 

oo 

and  was  well  approximated  by  the  following  simple  formula: 


6 = - (l  - |b  l/R  ) . 
4 o o 


Not  only  is  the  relationship  between  6 and  j B ^ | / R^  rather  insensitive 

to  geometry,  but  also  the  first-order  signal  statistics  are  not  very 

sensitive  to  that  relationship.  For  first-order  statistics,  the  main 

effect  of  diffraction  is  a simultaneous  rotation  (increasing  £>)  and 

circularization  (decreasing  |B  | /R  ) of  the  correlation  ellipse  illus- 

o o 

trated  in  Figure  1,  as  the  scattered  wave  propagates  with  increasing  Z in 

Eq.  (23)  . What  is  important  is  the  rate  at  which  this  approach  to  Rice 

statistics  (|b  |/R  =0  and  6 =45°)  takes  place  rather  than  the  precise 

o o 

relation  between  ellipse  elongation  and  eccentricity.  Thus,  we  have 

proceeded  on  the  basis  of  Eq.  (24) , which,  when  combined  with  Eq.  (7) , 

2 

yields  the  following  transcendental  equation  relating  |b^|/Rq  t0  °X^R0: 


lBo!  "A  j!o 

R 2 \ R 

o \ o 


= 1-  2r 


Eq.  (25)  is  fit  very  veil,  especially  in  the  rear  and  intermediate 
zones,  by  the  following  approximation: 


= 0.9811  - 1.55  — 
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Inserting  Eqs . (7)  and  (26)  into  Eq.  (5),  we  obtain  the  following  very 
useful  formula  for  scintillation  index: 


( / A / 2\  r / 2 

/ax\  2 /°x\  + R2  1 .98  - 1 . 52|  x 

S,  = 4R  — -4R  — o \ — 

4 o\R  / o\R  / \R 

\o/  \ o / \ o, 


While  containing  some  approximations  in  its  third  term,  Eq.  (27)  is  a 
considerable  improvement  over  the  formula  of  Briggs  and  Parkin  for 
[Eq.  (6)],  which  totally  neglected  both  the  second  and  third  terms.  The 
improvement  increases  as  the  scattering  becomes  more  intense  (increasing 


One  reason  for  separating  the  first  two  terms  in  Eq.  (27)  into  two 

factors  is  that  R is  amenable  to  a simple  multiple-scatter  correction, 
o 

The  theory  discussed  so  far  is  limited  to  the  regime  of  weak,  single 
scatter  for  the  entire  scattering  region  (i.e.,  Rq  « 1).  However,  an 
important  effect  of  multiple  scatter,  namely  extinction  of  the  non- 
scattered  signal  component  (represented  by  (E  in  Figure  1),  can  be 
accounted  for  simply  by  requiri.  g conservation  of  total  power  flux  while 
applying  the  first  Born  approximation  to  elemental  layers  within  the 
scattering  region.  From  the  work  of  Fejer  (1953)  , Bramley  (1955)  , 
Fremouw  (1968),  and  Rino  (1975),  it  follows  that  extinction  can  be 
accounted  for  by  replacing  the  single-scattering  coefficient  calculated 
from  Eq.  (16)  by  a multiple-scattering  coefficient  defined  by 


/ ~R  o 

R = 1 - e ° 
o 


It  is  R*  that  is  used  (in  place  of  R ) in  Eq.  (27)  for  calculation 
o ° 

rf  scintillation  index,  which  is  output  as  S4,  in  RFMOD ; this  multiple- 
scattering coefficient  (R  ) is  itself  output  fro  RFMOD  as  RZRO . An 
& o 
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important  task  for  the  second  contract  year  is  to  find  an  analogous 
multiple-scatter  correction  for  cj2/R  and  the  related  quantities  |b  |/R 

r X o o o 

and  6.  Presently,  the  latter  two  parameters  are  calculated  from  Eqs . 

(26)  and  (24),  respectively,  on  the  basis  of  Eq.  (23).  They  are  output 

from  RFMOD  as  BMAGR  and  DELTA,  respectively,  for  use  as  inputs  to  the 

second  channel-model  code,  DIST . A single-scatter  approximation  warning 

is  printed  if  R 5 1. 

o 

Once  R,  |b  | / R , and  6 are  calculated,  then  application  of  the 
o o o 

Gaussian-statistics  hypothesis  permits  calculation  of  any  desired  first- 
order,  signal-statistical  parameter.  The  calculation  of  the  intensity- 
scintillation  index, 

, * «»2  - uW>1/2 

4 (a2>  <29) 

as  described  above,  is  an  example.  An  analogous  phase-scintillation 
index  may  be  defined  as 


$ 


* <(*  . <$>)2>1/2 

rms 

= (<$2>  - <*>V/7' 


(30) 

(31) 


Given  the  phase  reference  chosen  in  Figure  1,  it  is  well  known  (e.g., 

Beckman  and  Spizzichino,  1963)  that  under  the  postulate  of  Gaussian 

statistics,  the  joint  probability  density  function  (pdf)  for  the  real, 

E , and  imaginary,  E , parts  of  the  received  signal's  complex  amplitude, 
x y 

E , is  given  by 

i r k -t,\2 


P (E  ,E  ) 
xy  x y 


2ttc  a 
x y 


2A/T77 


exp 


2(1  - P ) 


- 2p 


KE  - T])E 

x y 

o o 
x y 


+ 


(32) 
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9 / / 

where  p = o /a  0 and  where  T|  may  be  calculated  from  = V 1 * R an^ 
xy  x y ° 

the  o's  from  Eqs . (7),  (8),  and  (9).  The  individual  pdf's  for  amplitude 
and  for  phase  can  be  calculated  by  changing  variables  in  hq.  (32)  from 
component  to  polar  form  and  then  integrating  over  the  other  variab  e. 

That  is, 

2TT 

p ( A)  = \ A p (A  cos  4,  A sin  4)  d4  (33) 

Jo  xy 

CD 

and  p ( 4)  =1  A p (A  cos  4,  A sin  4)  dA  . (34) 

Jo  xy 

Analytic  solutions  to  the  integrals  in  Eqs.  (33)  and  (34)  have  been 

known  for  some  time  in  special  cases  (e.g.,  for  Rice  statistics,  in  which 

a =0  and  o = a ) , but  not  so  for  the  generalized  Gaus-'ian  statistics 
xy  x y 

needed  to  describe  scintillation.  Recently,  however,  Hatfield  (1975) 
has  found  the  following  general  analytic  solution  for  Eq.  (34): 


exp 


P(4)  = 


-v*> 

a 

_x 

1.2(1  - P2)  J I 


2TTCT  a A - P2  f 2(4) 
x y 1 


2 

- f2(0 

(1  - p ) exp 

2(1  - P2) 

. W1  - p2) 

+ f 2(4)  ^ erfc 


- f2(«) 
>4(1  - p2)j 


(35) 


where  f (4)  = \/cos2  4/o^  - 2P  cos  4 sin  4/(o^o^)  + sin 


2 2 
4/a 

y 


(36) 


2 

f (4)  = T]  [cos  4/a  - P sin  4/(a  a )j/f  (4) 

2 x x x y i 


and 


(37) 


U 'fortunately , there  is  no  known  analytic  solution  for  the  moments  of 
the  distribution  in  Eq.  (35).  Therefore,  the  distribution  itself  has 
been  incorporated  in  RFMOD,  together  with  a numerical  integration  routine 
for  calculating  its  first  and  second  moments,  as  follows: 


and 


2tt 

<*>  = ( 5 P (5)  dt 
2n 


p ($)  d$ 


(38) 


(39) 


The  results  of  Eqs . (38)  and  (39)  are  then  inserted  into  Eq.  (31)  for 
calculation  of  the  phase-scintillation  index,  which  is  output  as  PRMS . 

As  described  above,  the  signal-statistical  outputs  of  RFMOD  are 

r'  Ib  l/R  6 S and  $ . The  last  two  outputs  constitute  an  intensity- 

o’  1 o o’  ’ 4’  rms 

scintillation  index  and  a phase-scintillation  index,  respectively,  and 
represent  our  judgment  as  to  the  first-order  statistical  parameters  most 
likely  to  be  sought  by  a user  for  synoptic  purposes.  The  other  three 
outputs  are  provided  for  use  in  posing  more  detailed  questions  about 
first-order  signal  statistics  in  particular  cases,  by  means  of  the 
DIST  code 

The  DIST  code  is  essentially  a device  for  calculating  and  outputing 

the  amplitude  and  phase  pdfs  from  inputs  of  r\  IbJ/R^  and  6 • For 
phase,  it  is  simply  a matter  of  evaluating  Eq.  (35)  for  the  given  input 
parameters  and  outputing  the  result  as  P(PHI) . Unfortunately,  there  is 
no  known  analytical  solution  to  Eq.  (33)  analogous  to  Hatfield’s  phase 
distribution.  Thus,  it  is  necessary  to  perform  the  phase  integration 
over  the  joint  distribution  (32)  numerically.  This  integration  is  per- 
formed in  DIST,  and  the  result  is  output  as  P(A)  versus  A (the  latter  in 
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units  of  undisturbed  amplitude).  In  addition,  the  following  integration 


also  is  performed  numerically: 

V 

P(A  < V)  = j P ( A)  dA  . (40) 

t) 

This  resulting  cumulative  probability  function  (cpf)  of  amplitude  is 
then  output  as  CPF  versus  20  log  V (labeled  dB) . Since  the  independent 
variable  is  output  in  decibels,  CPF  may  be  thought  of  as  the  cumulative 
distribution  of  intensity  (I  = A ) as  well  as  of  amplitude  (A) , it 
represents  the  link  margin  needed  to  mitigate  scintillation  fading  to  a 
desired  reliability  threshold. 


Ill  SCINTILLATION  MORPHOLOGY  MODELED 


The  major  work  performed  on  this  contract  was  development  and 
implementation  of  the  generalized  scattering  theory  described  in  Section 
II.  The  other  important  aspect  of  modeling  the  trans ionospheric  scintil- 
lation channel  is  description  of  the  ionosphere's  scattering  structure  as 
a function  of  several  variables.  The  structural  characteristics  that 
must  be  described  pertain  to  the  strength,  size,  and  shape  of  the  radio- 
wave-scattering  irregularities  and  the  thickness  of  the  scattering  layer. 
All  the  relevant  parameters  appear  in  Eq.  (16).  The  shape  of  the  irregu- 
larities is  described  by  the  ionospheric  spatial  spectrum  ^Eqs.  (11)  and 
(11/ J,  which  is  characterized  mainly  by  the  power-law  spectral  index,  v, 
and  by  the  ratio,  a,  of  scale-size  parameter  along  the  geomagnetic  field 
to  that  across  it.  The  size  of  the  irregularities  is  scaled  by  the 
cross-field  scale-size  parameter,  a.  The  strength  of  the  irregularities 
is  described  by  the  spatial  variance  of  electron  density  ((AN)2),  and 
the  thickness  of  the  scattering  region  is  L. 

In-situ  measurements  of  the  ionospheric  spatial  spectrum  (Dyson  et 
al.,  1974)  show  a remarkably  consistent  spectral  index,  corresponding 
very  nearly  to  the  Kolmorov  value  of  4/3  for  v,  under  most  ionospheric 
conditions;  we  have  adopted  this  value  as  a constant  in  the  model.  The 
irregularity  axial  ratio,  a,  also  has  been  taken  as  constant,  with  a 
value  of  10,  as  a result  of  the  data  review  performed  by  Fremouw  and 
Rino  (1973).  There  is  no  direct  measurement  of  the  outer  scale,  a ; it  is 
simply  known  to  be  large  compared  with  a km  or  so  (Rufenach,  1972;  Dyson 
et  al.,  1974;  Sagalyn  et  al.,  1974).  We  have  found  that  a = 3 km  pro- 
duces a good  fit  to  the  amplitude-scintillation  data  analyzed  by  Rino, 
Livingston,  and  Whitney  (1975) , and  this  value  is  currently  employed  as 
a constant  in  the  model. 
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The  remaining  two  ionospheric  parameters,  ((AN)^  and  L,  appear  in 
the  scintillation  calculations  only  as  a product.  We  have  set  L = 100  km 
as  a representative  value  and  have  expressed  the  behavior  of  scintillation 
in  terms  of  the  following  empirical  model  for  the  rms  electron-density 
variation,  ((AN)2 


<(AN)2)1/2  = AN  (X,  T , D , R)  + AN  (X,T)  + AN  (X,T,K  ) 
e m n p 


+ AN  (X ,T,K  ) el/m 
a p 


where  X = geomagnetic  latitude  in  degrees 

T = time  of  day  in  hours, 

D = day  of  the  year  out  of  365, 

R = sunspot  number, 

K = planetary  magnetic  index. 

P 


(41) 


The  terms  in  Eq.  (41)  are 


9 

o 
»— 1 

+ 
t 

/t  + 1.5V 

= (1.2  X 10  )(1  + 0.05R) 

1 * °'25  coa  91.25  )\ 

exp 

[\  ‘ n 

+ exp 


T - 22.5 


-1 

/,  \2 

/ x \ 

exp 

"\12/ 

J 

. 

3 if  T £ 22.5 
TD  6 if  T > 22.5  (42) 


m 


r A \ 21 

/ rrT\ 

/X  - 32. 5\ 

^1  + 0.33  cos  — ) 

exn 

' 15 

L N ' - J 

(43) 


AN  = (2.0  X 10  ) 
h 


1 + erf 


68.5  + 1.82  Kp  + (5  + 0.55  Kp>  Cos  (ht/12) 


18.7  - 0.36  Kp  - (1  + 0.11  Kp)  cos(rrT/12) 


(44) 
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Eqs . (42)  and  (43)  describe,  respectively,  the  equatorial  and  mid- 
latitude behavior  of  scintillation,  and  the  high-latitude  behavior  is 
described  by  Eqs.  (44)  and  (45).  All  four  terms  of  Eq.  (41)  contain  the 
well-established  nighttime  maximum  of  scintillation.  The  equatorial  term 
also  includes  a seasonal  variation,  displaying  equinoctial  maxima.  The 
most  complicated  behavior  occurs  at  high  latitudes,  above  the  scintil- 
lation boundary"  (Aarons  et  al . , 1969).  Equation  (44)  essentially 
describes  the  behavior  of  the  boundary  location  as  a function  of  time  of 
day  and  geomagnetic  activity,  and  Eq.  (45)  accounts  for  scin' illation 
that  is  directly  related  to  auroral  activity. 

In  the  Fremouw-and-Rino  1973  model,  the  relationship  of  scintilla- 
tion to  geomagnetic  activity  was  not  included.  The  only  solar-terrestrial 
observable  employed  in  the  earlier  model  was  the  sunspot  number,  which 
was  intended  as  a measure  of  epoch  during  a solar  cycle.  However,  it  is 
clear  that  the  scintillation  boundary  migrates  equatorward,  along  with 
the  auroral  oval,  with  increasing  planetary  magnetic  activity  (Aarons 
et  »!.,  1963;  Evans,  1974).  This  was  known  at  the  time  the  1973  model 
was  developed,  but  the  dependence  was  excluded  to  avoid  double  accounting 
for  other  aspects  of  scintillation  behavior,  such  as  seasonal  and  solar- 
cycle  variation. 

Pope  (1974)  modified  our  high-latitude  term  to  include  a de- 
pendence without  checking  on  the  possibility  of  seasonal  double  accounting. 

We  have  had  five  years  of  K values  plotted  and  have  found  no  troublesome 

P 

seasonal  dependence.  We  have,  therefore,  adopted  Pope  s basic  revision 

of  replacing  the  sunspot-number  dependence  in  our  old  high-latitude  term 

with  a K dependence.  We  have  not  used  his  exact  formulation,  however, 

P 

which  deleted  our  earlier  description  of  the  variation  in  statistical 
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width  of  the  scintillation  boundary  with  its  location.  We  find  that  there 
is  such  a variable-width  behavior  in  the  in-situ  data  of  Sagalyn  et  al. 
(1974) ; it  is  described  by  the  denominator  of  the  error-function  argument 
in  Eq . (44) . 

Having  incorporated  a K dependence  in  Eq.  (44),  it  was  a simple 

extension  to  include  the  K -dependent  migration  of  the  auroral  oval  in 

P 

Eq.  (45)  . A more  accurate  index  to  use  might  be  the  auroral-electrojet 

index,  AE , but  we  have  elected  to  use  only  for  simplicity.  The  two 

magnetic  indices,  K and  AE , are  rather  highly  correlated,  as  shown  in 

P 

Figure  4.  To  avoid  the  necessity  for  a user  to  specify  two  magnetic 

indices  rather  than  just  one,  we  have  elected  to  exclude  AE  from  the 

model  unless  later  analysis  of  scintillation  data  shows  that  including 

it  would  improve  the  model's  accuracy  significantly.  In  a similar  vein, 

we  iiave  avoided  including  a geomagnetic  dependence  in  the  equatorial 

term,  retaining  instead  the  original  suns pot -number  dependence.  Further 

experimental  work  is  necessary  before  a reliable  geomagnetic  dependence 

can  be  incorporated  in  the  equatorial  term;  a likely  candidate  for  a 

magnetic  index  is  D , the  seasonal  and  solar-cycle  behaviors  of  which 

st 

are  similar  to  those  of  scintillation. 

In  addition  to  the  above  modifications  to  the  1973  model,  some  of  the 
constants  employed  have  been  changed.  These  changes  stem  partly  from 
improved  calibration  of  some  data  sets  used  in  development  of  the  original 
model.  The  data  sets  in  question  (Aarons  et  al.,  1963;  Koster,  1968; 

Joint  Satellite  Studies  Group,  1968)  were  presented  in  terms  of  a hand- 
scaled  scintillation  index,  SI  (Whitney  and  Malik,  1968),  which  only 
recently  (Whitney,  1974)  has  been  calibrated  reasonably  accurately  to  S^. 
It  is  still  to  be  noted  that  the  relation  between  SI  and  S^  must  be 
scintillation-rate  dependent,  and  data  characterized  directly  by  S^ 
are  strongly  preferred  for  modeling. 
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FIGURE  4 SCATTER  PLOT  SHOWING  DEGREE  OF  CORRELATION  BETWEEN  AURORAL  ELECTROJECT  INDEX 
AE,  AND  PLANETARY  MAGNETIC  INDEX,  Kp 


Figures  5 through  8 show  the  degree  of  fit  between  the  model 

now  in  RFMOD  and  several  data  sets  that  were  presented  by  the  observers 

either  in  terms  of  S or  in  terms  of  an  empirical  index  such  as  SI, 

4 

which  we  then  converted  to  C . Figures  9 and  10  compare  the  amplitude 

4 

pdf  and  cdf  calculated  by  means  of  DIST  with  corresponding  histograms 
developed  from  ATS  5 data  provided  by  Aarons  (private  communication, 
1974).  The  corresponding  pdf  for  phase  is  shown  in  Figure  11;  no  phase 
. data  exist  at  this  time  for  comparisot  . 

IV  THE  CHANNEL -MODEL  CODES 

Of  the  two  computer  programs,  RFMOD  and  DIST,  the  former  is  by 
far  the  more  complicated.  It  is  a modification  of  the  program  BPMOD, 
described  by  de  la  Beaujardiere  and  McNiel  (1971)  , and  consists  of  the 
routines  listed  in  Table  1.  Its  very  simplified  flow  chart  is  shown  in 
Figure  12.  Using  RFMOD  consists  essentially  of  preparing  input  cards, 
of  which  there  are  three  types:  (1)  a title  card,  (2)  namelist  cards, 

and  (3)  a changing-parameter(s)  card.  Three  sample  sets  of  input  cards 
are  shown  in  Figure  13,  and  each  type  of  card  is  described  below. 

The  title  card  is  used  to  assign  any  desired  title  of  up  to  80 
characters  to  a calculation  case  or  series  of  cases.  It  is  read  by 
the  driver  program,  RFMOD,  into  the  array  LBL  with  a 20A4  format,  and 
printed  in  the  output  by  the  subroutine  IMPRIME ; it  ha-  no  effect  what- 
soever on  the  rest  of  the  program.  Any  number  of  cases  can  be  run 
under  a common  title,  which  will  be  printed  in  the  output  for  each  case 

The  namelist  cards  are  used  to  input  the  fixed  parameters  for  a 
calculation  case.  They  are  read  by  subroutine  RF.ADIN,  by  means  of  tie 
namelist  INPAR.  The  user  gives  the  name  and  desired  value  for  each 
quantity  listed  above  the  dashed  line  in  Table  2 that  is  to  be  fixed. 
The  number  under  the  heading  K in  Table  2 is  an  index  used  in  the 
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FIGURE  9 COMPARISON  OF  CALCULATED  (smooth  curve) 
PROBABILITY  DENSITY  FUNCTION  FOR 
AMPLITUDE  WITH  OBSERVATIONS  (histogram) 
OF  ATS-5  TRANSMISSION  AT  137  MHz,  MADE 
BY  AARONS  (private  communication,  1974) 
THROUGH  THE  MID-LATITUDE  IONOSPHERE 
FROM  NEAR  BOSTON 


FIGURE  10  COMPARISON  OF  CALCULATED  (smooth  curve) 
CUMULATIVE  PROBABILITY  FUNCTION  (cpf) 
FOR  INTENSITY  (or  for  amplitude)  WITH 
OBSERVATIONS  (histogram)  FOR  CASE  IN 
FIGURE  9.  (Calculations  of  cpf  are  intended 
for  use  in  assessing  fade-margin  requirements.) 
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Table  1 


ROUTINES  IN  PROGRAM  RFMOD 


Name 

Function 

BLOCK  DATA 

Initializes  common  parameters. 

SUBROUTINE  READIN, 
with  ENTRY  INC 

Performs  program  inputting,  and  increments  the 
changing  parameter(s) . 

SUBROUTINE  SCINT 

Controls  the  basic  scintillation  calculation  and 
provides  quantities  needed  for  SUBROUTINE  SIGS. 

SUBROUTINE  SIGS 

Performs  the  bpsic  scintillation  calculations. 

FUNCTION  RMSDN 

Contains  the  worldwide  model  for  strength  of 
F-layer  irregularities  used  in  the  scintillation 
calculation . 

FUNCTION  RFBPT 

Evaluations  the  diffraction  integral  by  an 
algebraic  approximation. 

SUBROUTINE  FINDZ 

Computes  the  reduced  distance  from  scattering 
layer  to  receiver  for  use  in  evaluating  Fresnel- 
zone  radius. 

SUBROUTINE  MGFLD 

Contains  an  earth-centered,  axially  tapped, 
dipole  model  of  the  geomagnetic  field  and  cal- 
culates the  geomagnetic  coordinates  of,  and  the 
field  components  at,  a specified  point. 

SUBROUTINE  AZNSIDE 

Given  the  latitude  and  longitude  of  two  points, 
A and  B,  computes  the  azimuth  of  B as  seen  from 
A and  the  great-circle  arc  AB. 

SUBROUTINE  COORD 

Calculates  coordinates  of  a desired  point  on  a 
specified  great-circle  arc. 

FUNCTION  ERF 

Computes  the  Error  Function. 

SUBROUTINE  CPRMS 

Controls  calculation  of  rms  phase. 

FUNCTION  PHAS, 
with  ENTRY  PRT 

Evaluates  the  Hatfield  phase  distribution  of 
a scintillating  signal. 

FUNCTION  FI,  with 
ENTRY  F2 

Provides  integrand  functions  for  evaluating 
moments  of  phase  distribution. 

SUBROUTINE  SIMPA 

Adaptively  performs  numerical  integration 
using  Simpson's  rule. 

SUBROUTINE  IMPRIME 

Performs  program  outputting. 
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DISPLAY  THE  VARIOUS  PARAMETERS  FOR  THIS 
CASE,  AND  THE  VALUES  OF  THE  SCINTILLATION  INDICES 
AND  THE  PROBABILITY-ELLIPSE  PARAMETERS 


LA-1079-28R 


FIGURE  12  SIMPLIFIED  FLOW  CHART  OF  PROGRAM  RFMOD 
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Table  2 

DATA-CARD  VARIABLES 


- 

K Variable  Name 

Unit 

Definition 

1 FREO 

Megahertz 

Frequency 

2 SSN 

None  (floating  point) 

Sunspot  number 

3 DAY 

None  (floating  point) 

Day  of  year 

4 TIME 

Decimal  hours 

Meridian  time  at  iono- 

spheric  penetration  point 

5 RLAT 

Decimal  degrees 

Receiver  latitude 

6 RLON 

Decimal  degrees 

Receiver  longitude 

7 HR 

Me  te  rs 

Receiver  altitude 

8 TLAT 

Decimal  degrees 

Transmitter  latitude 

9 TLON 

Decimal  degrees 

Transmitter  longitude 

10  HT 

Me  ters 

Transmitter  altitude 

11  FKP 

Decimal  value  0.0 

Planetary  magnetic 

through  9.0 

index,  Kp 

11 

None  (value  1 or  0) 

Data  flag  (see  text) 

12  RCRD 

Decimal  degrees 

Receiver  coordinates 

13  TCRD 

Decimal  degrees 

Transmitter  coordinates 

common  VARPAR  to 

identify  the  parameters. 

The  user  may  list  the  fixed 

parameter?  in  any  order  on  the  namelist  cards. 

A parameter  input  on  a namelist  card  will  keep  the  specified  value 
until  changed  by  means  of  a new  namelist  card  or  until  input  on  a changing- 
parameter(s)  card.  The  first  column  of  each  namelist  card  must  be  blank, 
and  the  first  characters  on  the  first  namelist  card  must  be  &INPAR . The 
last  namelist  card  in  a sequence  must  terminate  with  the  characters  &END. 


The  changing-parameter(s)  card  is  employed  by  the  user  to  select  the 
independent  variable(s)  for  a given  calculation  case.  He  may  choose 
either  of  two  program  modes.  In  Mode  1,  he  chooses  any  one  of  the  param- 
eters listed  in  Table  2 for  which  K is  between  1 and  11,  inclusive,  and 
inputs  its  name  on  the  changing-parameter  card.  In  this  mode,  the  first- 
order  signal-statistical  parameters,  including  and  ^rnig>  are  calculated 
as  functions  of  the  single  changing  parameter.  In  Figure  13,  the  cards 
identified  as  Set  1 are  for  two  sequential  cases  run  in  Mode  1,  under  a 
common  tLtle. 

Mode  2 permits  the  user  to  increment  the  position  of  either  the  re- 
ceiver or  the  transmitter  along  an  arc  of  a great  circle.  In  this  mode, 
both  the  latitude  and  longitude  of  the  receiver  or  of  the  transmitter  be- 
come changing  parameters,  and  the  signal-statistical  parameters  are  cal- 
culated for  equidistant  points  along  the  arc.  Mode  2 is  selected  by  namirg 
RCRD  or  TCRD  on  the  changing-parameters  card  (see  Table  2).  In  Figure  13, 
the  cards  identified  as  Set  2 are  for  a run  in  this  mode. 

The  required  contents  of  the  changing-parame ter( s)  card  are  listed 
in  Table  3.  Note  that  in  Mode  1,  if  the  final  value  of  a changing  param- 
eter is  smaller  than  the  initial  value,  the  increment  value  must  be  speci- 
fied as  negative.  Absence  of  the  required  negative  sign  in  this  situation 
will  cause  the  run  to  terminate,  and  an  error  message  will  be  printed. 

After  reading  a set  of  cards,  the  program  will  keep  the  same  title 
and  read  only  namelist  and  c nging-parame ter( s)  cards,  so  long  as  II  = 1. 
After  encountering  II  = 0,  it  will  read  a new  title  card.  The  program 
will  terminate  if  the  name  on  the  changing-parame ter ( s)  card  does  not 
correspond  to  one  of  those  numbered  1 through  13  in  Ta^e  2.  Thus,  to 
terminate  a run,  the  user  may  prepare  a changing-parameter  card  contain- 
ing an  unfamiliar  word,  such  as  the  FINI  appearing  in  Figure  13. 
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Table  3 


CONTENTS  OF  CHANGiNG-PARAMETER( S)  CARD 


Columns 

Mode-1  Content 

Mode-2  Content 

Symbol 
in  Program 

1 through  4 

Name 

Name 

Name 

9 through  16 

Initial  value 

Initial  latitude 

A1 

17  through  24 

Final  value 

Initial  longitude 

42 

25  through  32 

Increment  value 

Final  latitude 

A3 

33  through  40 

Final  longitude 

A4 

41  through  42 

-- 

Number  of  increments 

NA5 

For  model  development  and  testing,  we  employed  a plotting  routine 
in  conjunction  with  RFMOD  to  produce  graphs  such  as  those  in  Figures  5 
through  8.  The  available  vers  ion  of  RFMOD  provides  an  array  that  the  user 
may  employ  with  his  own  computer  peripherals.  In  the  absence  of  user 
modification,  the  array  will  be  output  in  tabular  form  along  with  identify- 
ing information,  as  shown  for  a test  case  in  Figure  14.  In  the  case  shown, 
the  changing  parameter  was  the  latitude  of  the  transmitter,  which  was 
taken  to  be  in  a nolar-orbiting  satellite  passing  over  a receiving  station 
near  Boston  (namely,  the  Sagamore  Hill  site  of  the  AFCRL  group  headed  by 
Aarons).  The  changing  parameter  is  followed  in  each  row  by  the  azimuth 
and  elevation  of  the  satellite,  the  three  probability-ellipse  parameters, 
and  then  by  the  calculated  scintillation  indices.  The  input  cards  for 
this  run  are  identified  as  Set  3 in  Figure  13. 

It  is  envisioned  that  a user  will  inspect  the  values  of  S4  and/or 
PRMS  predicted  for  his  communication  link  by  RFMOD  and  then  decide  if  he 
wants  more  detailed' signal-s tatistical  information  for  a particular 
situation.  In  the  example  of  Figure  14,  for  instance,  he  might  decide 


-LAYER  produced  vhf-uhe  scintillation 


SINGLfc-SCaTTEP  aPPPO*.  warning 

FIGURE  14  EXAMPLE  OF  OUTPUT  FROM  AVAILABL  VERSION  OF  PROGRAM  RFMOD 


thathewants  the  full  first-order  distributions  of  intensity  and  phase 
for  the  case  in  which  S4  on  his  channel  is  .538  and  PRMS  is  32.835  degrees. 
To  obtain  the  distributions  he  would  employ  DIST,  inputting  RZRO  = .371, 
BMAGR  = .924,  and  DELTA  = 3.419  degrees,  on  a single  card  in  F10.3  format. 
The  DIST  code  consists  of  the  routines  listed  in  Table  4. 


Table  4 

ROUTINES  IN  PROGRAM  DIST 


Name 

Function 

SUBROUTINE  PDFI 

Controls  calculation  and  output  of  probability 
density  function  for  amplitude,  by  means  of  one- 
dimensional integration  over  two-dimensional 
Gaussian  probability  density  function. 

SUBROUTINE  PARRAY 

Initializes  amplitude  array  for  SUBROUTINE  PDFI. 

SUBROUTINE  GAUP 

Controls  calculation  of  Gaussian  probability 
density  function. 

! FUNCTION  PR 

Computes  two-dimensional  Gaussian  probability 
density  function. 

FUNCTION  XINTGR 

Performs  integration  by  Simpson's  rule. 

SUBROUTINE  PDFP 

Controls  computation  and  output  of  probability 
density  function  for  phase,  calling  SUBROUTINE 
PHAS . 

SUBROUTINE  PHAS 

Evaluates  Hatfield  phase  distribution. 

FUNCTION  ERF 

Computes  the  Error  Function. 

SUBROUTINE  CPF I 

Calculates  and  outputs  the  cumulative  probability 
density  function  for  amplitude  or  intensity 
(expressed  in  dB)  by  performing  Simpson's  rule 
integration  over  probability  density  function 
for  amplitude. 
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Again,  for  modeling  we  employ  a plotting  routine  in  conjunction  with 
DIST  to  produce  graphs  such  as  those  in  9,  10,  and  11.  The  available 
version  of  DIST  simply  outputs  the  array  from  which  such  graphs  are  made, 
in  the  line-printer  formats  shown  in  Figures  15,  16,  and  17.  As  with 
RFMOD,  a user  may  wish  to  employ  DIST  in  conjunction  with  his  own  plotting 
routine  or  in  some  other  output  form.  Alternatively,  he  may  wish  to 
reformat  the  array  listing,  perhaps  outputting  fewer  values  and  possibly 
presenting  CPF  in  terms  of  percent  or  decimal  link  reliability  (the  one's 
complement  of  CPF)  rather  than  directly  as  in  Figure  16. 

V CONCLUSION 

It  has  been  our  attempt  to  develop  a channel  model  that  may  be  used 
conveniently  for  calculating  first-order  signal  statistics,  expressed  in 
a format  relevant  to  systems  application.  The  resulting  codes,  RFMOD 
and  DIST,  represent  the  first-year  output  of  a broader  effort  to 
definitively  model  the  transionospheric  communication  channel.  These 
codes  are  to  be  augmented  by  routines  to  permit  calculation  of  second- 
order  statistics,  such  as  fluctuation  spectra  of  amplitude  and  phase,  as 
well  as  signal  correlation  at  spaced  receivers  and  across  bands  of 
frequency.  Moreover,  the  existing  version  of  RFMOD  is  likely  to  evolve 
over  the  next  year  in  two  important  ways. 

First,  it  is  intended  to  modify  SUBROUTINE  SIGS  in  RFMOD  to  fully 
account  for  multiple  scatter.  Theoretical  work  now  underway  is  hinting 
at  a possible  form  for  the  modification.  The  single-scatter  theory 
described  in  Section  II  of  this  report  describes  the  simultaneous  rota- 
tion (increasing  6)  and  circularization  (decreasing  |B0|/R0)  of  the 
first-order  probability  ellipse  for  the  received  signal's  complex  ampli- 
tude (illustrated  in  Figure  1),  as  the  wave  propagates  away  from  the 
scattered  medium.  In  the  single-scatter  limit,  the  ellipse  emerges  from 
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EXAMPLE  OF  AMPLITUDE  pdf  OUTPUT  FROM  AVAILABLE  VERSION  OF  PROGRAM  DIST 
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the  scattering  layer  highly  elongated  (|B0|/R0  « 1)  and  oriented  nearly 
vertically  (6  « 0)  on  the  complex  plane. 

Multiple-scatter  calculations  for  the  special  case  of  normal  inci- 
dence on  the  scattering  layer  show  that  some  circularization  (|B0|/R0 
rj  e“Ro),  but  no  rotation  (6  r;  0),  is  produced  Dy  the  rescat-erings.  This 
is  an  intuitively  understandable  result  and,  if  corroborated  by  similar 
calculations  for  off-normal  incidence,  is  expected  to  lead  to  a multiple- 
scatter  correction  for  |Bq|/R0  analogous  to  Eq.  (28)  for  RQ.  It  will 
presumably  also  lead  tc  a modification  of  Eq.  (24)  for  relating  the 
ellipse  orientation  to  its  elongation  at  the  receiver.  The  result  is 
certain  to  be  a reduction  in  the  third  term  of  Eq.  (5)  for  intensity- 
scintillation  index  as  the  fraction  of  flux  scattered  (Rq)  increases. 

It  will  presumably  limit  values  of  to  something  near  unity  and  improve 
agreement  between  calculated  and  observed  values  for  moderate  to  strong 
sc intillation. 

The  other  important  area  in  which  RFMOD  is  expected  to  evolve  is 

that  of  scintillation  morphology,  which  will  entail  updating  of  FUNCTION 

RMSDN.  As  more  definitive  scintillation  data  become  available,  the  model 

of  electron-density  irregularity  strength  contained  in  RMSDN  presumably 

will  be  improved  by  a process  that  may  be  viewed  as  calibration  [i.e., 

the  main  changes  should  be  adjustment  of  the  coefficients  in  Eqs.  (4^) 

through  (45)].  This  will  be  particularly  so  in  regions  where  the  model 

is  currently  calibrated  in  terms  of  ad  hoc  indices  such  as  SI  instead  of 

in  terms  of  S, . One  such  region  is  the  important  one  centered  on  the 
4 

geomagnetic  equator.  Useful  equatorial  data  are  expected  from  observa- 
tions of  the  NOAA  ionospheric  beacon  on  board  NASA's  ATS-6,  particularly 
measurements  of  S.  at  140  and  360  MHz.  During  periods  of  weak  to 
moderate  scintillation,  estimates  should  also  be  obtained  of 
at  140  MHz.  The  ATS-6  observations  also  will  provide  histograms 


of  amplitude  (under  all  conditions))  and  of  phase  (under  conditions  of 
weak  to  moderate  scintillation)  for  testing  the  predictions  made  by  the 
DIST  code. 

Pending  the  improvements  described  above,  we  offer  RFMOD  and  DIST 
as  the  best  way  to  estimate  the  first-order  signal  statistics  to  be 
expected  at  the  output  of  a user's  transionospheric  communication  link, 
short  of  performing  a series  of  observations  over  that  particular  link. 

It  must  be  remembered  that  the  estimated  values  are  those  expected  for 
average  ionospheric  conditions  in  the  observing  situation  (frequency, 
transmitter  and  receiver  location,  time  of  day,  K^,  etc.)  specified  by 
the  user.  There  is  presumably  some  statistical  distribution  of  each 
parameter  estimated  about  its  calculated  value,  and  the  present  model 
does  not  account  for  such  statistical  variation.  Possibly  of  equal 
importance,  the  model  has  yet  to  be  subjected  to  testing  against  observa- 
tions above  400  MHz  (e.g.,  in  the  SHF  band). 


REFERENCES 


Aarons,  J.,  J.  P.  Mullen,  and  S.  Basu,  "Geomagnetic  Control  of  Satellite 
Scintillations,"  J.  Geophys.  Res. , Vol.  68,  No.  10,  pp.  3159-3168 
(1963). 

Aarons,  J.,  J.  P.  Mullen,  and  H.  E.  Whitney,  "The  Scintillation  Boundary," 
J.  Geophys.  Res. , Vol.  74,  No.  3,  pp.  884-889  (1969). 

Beckman,  P. , and  A.  Spizzichino,  The  Scattering  of  Electromagnetic  Waves 

from  Rough  Surfaces,  pp.  119-136  (Pergamon  Press,  New  York,  New  York, 
1963). 

Bramley,  E.  N. , "Some  Aspects  of  the  Rapid  Directional  Fluctuations  of 

Short  Radio  Waves  Reflected  at  the  Ionosphere,"  Proc.  IEEE,  Vol.  102, 
pp.  533-540  (1955). 

Briggs,  B.  H.,  and  I.  A.  Parkin,  "On  the  Variation  of  Radio  Star  and 

Satellite  Scintillations  with  Zenith  Angle,"  J.  Atmos.  Terr.  Phys. , 

VoJ.  25,  pp.  339-365  (1963). 

Budden,  K.  G. , "The  Amplitude  Fluctuations  of  the  Radio  Wave  Scattered 

from  a Thick  Ionospheric  Layer  with  Weak  Irregularities,"  J.  Atmos. 
Terr.  Phys. , Vol.  27,  pp.  155-172  (1965). 

de  la  Beaujardiere,  0.,  and  D.  McNiel,  "A  Fortran  Program  for  Calculating 
F- Layer-Produced  Scintillation,"  Final  Report,  SRI  Project  10 ’9, 
Contract  NAS5-21551,  Stanford  Research  Institute,  Menlo  Park, 
California  (December  1971). 

Dyson,  P.  L. , J.  P.  McClure,  and  W.  B.  Hanson,  "In  Situ  Measurements  of 

the  Spectral  Characteristics  of  F Region  Ionospheric  Irregularities," 
J.  Geophys.  Res..  Vol.  79,  No.  10,  pp.  1497-1502  (1  April  1974). 

Evans,  J.  V.,  Ed.,  "Millstone  Hill  Radar  Propagation  Study:  Scientific 

Results--Part  II,"  Technical  Report  509,  Lincoln  Laboratory, 
Massachusetts  Institute  of  Technology,  Lexington,  Massachusetts 
(12  November  1973). 


Fejer,  J,  A,  , "The  Diffraction  of  Waves  in  Passing  Through  an  Irregular 
Refracting  Medium,"  Proc.  Roy.  Soc.  A. , Vol.  220  (1953). 

Fremouw,  E,  J. , Ed.,  "An  Experiment  to  Measure  Trans-Atmospheric  Wideband 
Radio  Propagation  Parameters,"  Final  Report--Phase  I,  Contract 
DASA01-68-C-0104,  DASA  Report  DASA-2362-1,  SRI  Project  7161, 

Stanford  Research  Institute,  Menlo  Park,  California  (1969). 

Fremouw,  E.  J. , "Measurement  of  Phase  Variance  and  Autocorrelation  During 
Radio-Star  Visibility  Fades,"  J.  Geophys.  Res. , Vol.  73,  No.  11, 
pp.  3557-3564  (1968). 

Fremouw,  E.  J. , and  C.  L.  Rino,  "An  Empirical  Model  for  Average  F-Layer 

Scintillation  at  VHF/UHF,"  Radio  Science,  Vol.  8,  No.  3,  pp.  213-222 
(March  1973). 

Hatfield,  V.  E.,  "Non-Rician  Statistics  and  Their  Implications  for  Model- 
ing Effects  of  Scintillation  on  Communication  Channels,"  paper 
presented  at  NRL  Symposium  on  Effects  of  the  Ionosphere  on  Space 
Systems  and  Communication,  Washington,  D.C.,  20-22  January  1975. 

JSSG  [Joint  Satellite  Studies  Group],  "On  the  Latitude  Variation  of 

Scintillations  of  Ionosphere  Origin  in  Satellite  Signals,"  Planet . 
Space  Sci.  , Vol.  16,  No.  6,  pp.  775-781  (1968). 

Foster,  J.  R. , "Equatorial  Studies  of  the  VHF  Signal  Radiated  by 

Intelstat  II,  F-3,  1,  Ionospheric  Scintillation,"  Progress  Report  3, 

Contract  F61052-67-C00026,  University  of  Ghana-Legon,  Accra,  Ghana. 

Nakagami,  M. , "The  m-Distribution--A  General  Formula  of  Intensity 

Distribution  of  Rapid  Fading,"  Stat.  Methods  in  Radio  Wave  Prop., 

W.  C.  Hoffman,  ed.,  pp.  3-36  (Pergamon  Press,  Oxford,  1960). 

Pope,  J.  H.,  "High-Latitude  Ionospheric  Irregularity  Model,"  Radio 
Science.  Vol.  9,  No.  7,  pp.  576-682  (1974)  . 

Rj.no,  C.  L.  , "Analysis  of  Scintillation  Effects  on  Communication  and 
Radar  Systems,"  Bimonthly  Progress  Report  4,  SRI  Project  3399, 
Contract  DNA001- 74-C-0255,  Stanford  Research  Institute,  Menlo 
Park,  California  (15  January  1975). 

Rino,  C.  L.  , and  E..  J.  Fremouw,  "Statistics  for  Ionospherically  Diffracted 
VHF/UHF  Signals,"  Radio  Science,  Vol.  8,  No.  3,  pp.  223-233 
(March  1973a). 


48 


Rino,  C.  L. , and  E.  J.  Fremouw,  "Ionospheric  Scintillation  Studies," 

’ Final  Report,  Contract  NAS5-21891,  SRI  Project  2273,  Stanford 
Research  Institute,  Menlo  Park,  California  (November  1973b). 

Rino,  C.  L. , R.  C.  Livingston,  and  H.  E.  Whitney,  "Some  New  Results 
on  the  Statistics  of  Radio-Wave  Scintillation:  A.  Empirical 

Evidence  for  Gaussian  Statistics,"  accepted  for  publication  in 
J.  Geophys.  Res. , 1975. 

Rufenach,  C.  L. , "Power-Law  Wavenumber  Spectrum  Deduced  from  Ionospheric 
Scintillation  Observations,"  J.  Geophys.  Res.,  Vol.  77,  No.  25, 
pp.  4761-4772  (1  September  1972). 

Sagalyn,  R. , M.  Smiddy,  and  M.  Ahmed,  "High-Latitude  Irregularities  in 
the  Top  Side  Ionosphere  Based  on  ISIS  I Thermal  Ion  Probe  Data," 
j.  Geophys.  Res.,  Vol.  79,  No.  28,  pp.  4252-4261  (October  1974). 

Whitney,  H.  E.  , "Notes  on  the  Relationship  of  Scintillation  Index  to 
Probability  Distributions  and  Their  Uses  for  System  Design," 
Environmental  Research  Papers  No.  461,  Report  No.  AFCRL-TR-74-0004, 
Ionospheric  Physics  Laboratory,  Air  Force  Cambridge  Research 
Laboratories,  L.  G.  Hanscom  F'.eld,  Bedford,  Massachusetts  (January 
1974). 

Whitney,  H.  E. . and  C.  Malik,  "A  Proposed  Index  for  Measuring  Ionospheric 
Scintillation,"  Environmental  Research  Papers  No.  284,  Report  No. 
AFCRL-68-0138,  Ionospheric  Physics  Laboratory  Air  Force  Cambridge 
Research  Laboratories,  L.  G.  Hanscom  Field,  Bedford,  Massachusetts 
(March  1968). 


49 


